A Terradynamics of Legged Locomotion on Granular Media 
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The theories of aero- and hydrodynamics pre- 
dict animal movement and device design in air 
and water through the computation of lift, drag, 
and thrust forces. Although models of terrestrial 
legged locomotion have focused on interactions 
with solid ground, many animals move on sub- 
strates that flow in response to intrusion. How- 
ever, locomotor-ground interaction models on 
such flowable ground are often unavailable. We 
developed a force model for arbitrarily-shaped 
legs and bodies moving freely in granular media, 
and used this "terradynamics" to predict a small 
legged robot's locomotion on granular media us- 
ing various leg shapes and stride frequencies. Our 
study reveals a complex but generic dependence 
of stresses in granular media on intruder depth, 
orientation, and movement direction and gives in- 
sight into the effects of leg morphology and kine- 
matics on movement. 

The locomotion of animals (1) and devices (2-4) 
emerges from the effective interaction of bodies and/or 
appendages with an environment. For flying in air and 
swimming in water, there is a history of theoretical pre- 
dictive models (3-5) to describe the complex interactions 
between the locomotor and the surrounding fluids, based 
on the fundamental equations for fluid flow, the Navier- 
Stokes equations. These models have not only allowed 
understanding of the movement of a variety of aerial and 
aquatic organisms (5) [such as bacteria and spermatazoa 
(6), insects (7), birds (8), and fish and whales (9)] and 
their functional morphology, evolution, and ecology (9, 
10), but also advanced the engineering design of aircraft 
(3), marine vehicles (4), and flying (11) and swimming 
(12) robots. For running and walking on ground, studies 
using solid ground such as running tracks and treadmills 
have inspired general models (13, 14); building on these 
models, researchers have begun to apply dynamical sys- 
tems theory (15). In these studies, the leg-ground inter- 
action was often approximated as a point contact on a 
rigid, flat, and nonslip ground (13-15). 

Many small legged animals (16-19) [and increasingly 
robots (20-23)] face the challenges of moving on natural 
substrates such as sand (16, 17, 21), gravel (16, 20), rub- 
ble (20), soil (20, 22),mud (17, 20), snow (18, 20), grass 
(20, 22), and leaf litter (19, 20, 22), which, unlike solid 
ground, can flow during movement when a yield stress is 
exceeded. The complexity of the interactions with such 
flowable ground may rival or even exceed that during 



movement in fluids. For example, recent studies of legged 
animals (16) and robots (21) moving on granular media 
[collections of particles (24)] such as sand and gravel (Fig. 
1, A and B) have demonstrated that at an instant of time 
during a step, each element of a leg moves through the 
substrate at a specific depth, orientation, and movement 
direction, all of which can change over time (16, 21). Fur- 
thermore, the leg interacts with a material that can dis- 
play both solid-like and fluid-like features (24) (Fig. 1, C 
and D). Compared to the theories of aero- and hydrody- 
namics, predictive models are less well developed for cal- 
culating forces and predicting legged locomotion on such 
flowable ground (16, 21). Research in the field of terrame- 
chanics (2) has advanced the mobility of off-road vehicles 
on flowable ground such as sand and soil. These models 
were developed for large wheeled and tracked vehicles, 
which sink only slightly into the substrate (2, 25). Thus, 
in terramechanical models, interaction with the ground 
is approximated as the indentation of a horizontal, flat, 
rectangular plate (2, 25). It was a breakdown of this flat- 
plate approximation, however, that led to overpredicted 




FIG. 1. Examples of legged locomotion on flowable ground. 
(A) A zebra-tailed lizard running on sand (16). (B) A bi- 
ologically inspired RHex robot (22) walking on dirt [Photo 
credit: Galen Clark Haynes, Aaron M. Johnson, and Daniel 
E. Koditschek, University of Pennsylvania]. Dashed boxes 
in (A) and (B) indicate the regions of leg-ground interaction 
shown in (C) and (D). Schematic of leg-ground interaction for 
(C) a hind foot of a zebra-tailed lizard (16) and (D) a c-leg of 
a RHex robot (21) during a step on granular media. Dashed, 
solid, and dotted tracings are leg positions at early, mid-, and 
late stance. Bars and arrows indicate local orientations and 
movement directions of leg elements. The gray area is the 
granular substrate. 
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speeds for small vehicles such as the Mars rovers, whose 
small wheels have substantially curved ground contact 
interfaces (25). Because leg-ground interaction on flow- 
able substrates is a more diverse, complex, and dynamic 
process (16, 21) than the flat-plate indentation, terrame- 
chanics is not expected to apply to legged locomotors on 
flowable ground (2). 

Granular materials such as sand and gravel (24), home 
to a variety of small desert animals (16, 17, 26), have 
proved to be a promising model substrate for studying 
legged locomotion on flowable ground (16, 21). Despite 
their diversity in particle size, shape, density, friction, 
polydispersity, and compaction (24), dry granular me- 
dia are relatively simple as compared to media such as 
soil and mud, because granular particles interact purely 
through dissipative, repulsive contact forces and have no 
cohesion (24). In addition, the penetration resistance of 
granular media can be repeatably controlled using lab- 
oratory apparatus such as an air-fluidized bed (16, 21, 
26). 

To begin to create a "terradynamics" that allows the 
prediction of legged locomotion on a flowable ground, we 
hypothesized that the net forces on a leg (or a body) 
moving in a granular medium in the vertical plane could 
be approximated by the linear superposition of resistive 
forces on infinitesimal leg (or body) elements. Our hy- 
pothesis was inspired by our recent success in applying 
the methods of resistive force theory (6) to predict the 
forces and movement during the limbless locomotion of a 
lizard swimming in sand (26) and in describing the drag 
(26, 27) and lift (27) on simple objects moving in gran- 
ular media at fixed depths. In these studies, the linear 
superposition was valid for intruders moving in granu- 
lar media in the horizontal plane at low enough speeds 
[for example, < 0.5 m/s for 0.3-mm glass particles (26)], 
where intrusion forces are dominated by particle friction 
(insensitive to speed) and non-inertial (26). However, it 
was unclear whether linear superposition could apply to 
legs (or bodies) of complex morphology and kinematics 
moving in the vertical plane. 

To measure resistive forces for leg elements, we moved 
a thin rigid plate (of area A) in granular media in the 
vertical plane at 1 cm/s and measured lift f z and drag 
f x on the plate (in the continuously yielding regime). We 
determined vertical and horizontal stresses a ZjX = f ZjX /A 
as a function of the plate's depth \z\ below the surface, 
angle of attack /3, and angle of intrusion 7 (Fig. 2 A and 
movie SI) (28). To test the generality of our resistive 
force model, we used three dry granular media of vari- 
ous particle size, shape, density, and friction, prepared 
into flat, naturally occurring, loosely and closely packed 
states (16, 21, 26) (supplementary text section 2, fig. S3, 
and table SI). Slightly polydispersed near-spherical glass 
particles 0.3 and 3 mm in diameter [covering the particle 
size range of natural dry sand (~0.1 to ~1 mm) (29)] and 
rounded, slightly kidney-shaped poppy seeds (0.7 mm in 



diameter) allowed us to probe general principles for nat- 
urally occurring granular media of high sphericity and 
roundness [such as Ottawa sand (30)]. We discuss at the 
end of the paper possible effects of particle nonsphericity 
and angularity also found in many natural sands (30). 

In all media tested, we observed that for all attack 
angles f3 and intrusion angles 7, a ZyX were nearly propor- 
tional to depth \z\ when the plate was fully submerged 
and far from the bottom of the container (Fig. 2B and 
movie SI). This is because friction-dominated forces are 
proportional to the hydrostatic-like pressure in granu- 




FIG. 2. Measurement of resistive forces in granular media 
in the vertical plane using a plate element (28). (A) Lift f z 
(blue arrow) and drag f x (green arrow) on a thin rigid plate 
of area A moving in granular media (gray area) in the verti- 
cal plane at speed v = 1 cm/s were measured as a function 
of the plate's depth \z\ below the surface, angle of attack /3, 
and angle of intrusion 7. The granular media were fluidized 
(and then compacted when a closely packed state was pre- 
pared) using an air-fluidized bed (26) before each intrusion 
(7 > 0) and extraction (7 < 0). g is gravitational accelera- 
tion. (B) Vertical (blue curve) and horizontal (green curve) 
stresses a z , x = f z ,x/A versus \z\ for representative intrusion 
and extraction using (^,7) = ±(7r/6, 7r/4) (movie SI). Blue 
and green dashed lines are linear fits with zero intercept over 
intermediate depths at which the plate was fully submerged 
and far from the bottom of the container. Horizontal ar- 
rows on top indicate the range of leg depths in Figs. 3 and 
4. (C) Vertical and (D) horizontal stresses per unit depth 
a z , x [slopes of dashed fit lines in (B)] versus f3 and 7. Plate 
schematics with arrows denote representative orientations and 
movement directions. Black curves indicate where a z , x = 0. 
The shaded areas indicate where a z (or a x ) is not opposing 
the plate's vertical (or horizontal) velocity. Circles indicate 
a z , x from data shown in (B). Arrows above and below the 
color bars indicate directions of a z , x for positive and negative 
values. 
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lar media. Therefore, we modeled the hydrostatic-like 

stresses as 
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where a z ^ x are vertical and horizontal stresses per unit 
depth (slopes of dashed fit lines in Fig. 2B). We found 
that in all media tested, a ZjX depended sensitively on 
both attack angle (3 and intrusion angle 7 (Fig. 2, C and 
D, fig. S4, and additional data table S5). a z (or a x ) 
was opposing the plate's vertical (or horizontal) velocity 
for most but, counterintuitively, not all (3 and 7 (excep- 
tions are indicated by the shaded regions). For almost all 
attack angles /?, a ZjX had larger magnitudes (|a^ >£C |) for 
intrusion angle 7 > than for 7 < 0; i.e., it was harder 
to push the plate into granular media than to extract it. 
For all intrusion angles 7 except 7 = ±7r/2, \a ZjX \ were 
asymmetric to attack angles (3 = and (3 = ±7r/2; i.e., 
only when the plate moved vertically were stress mag- 
nitudes the same for vertically or horizontally mirrored 
orientations (e.g., (3 = ±7r/6). These asymmetries are a 
result of gravity breaking symmetry in the vertical plane 
and differ from the case in the horizontal plane (26). Our 
resistive force measurements are an advance from previ- 
ous force models based on the flat-plate approximation 
used in many terramechanical models (2, 25), which cap- 
ture only the dependence of stresses on intruder depth, 
but not on its orientation or movement direction (sup- 
plementary text section 1 and fig. SI). Despite their 
different magnitudes and subtle differences in shape, the 
overall profiles of stresses (per unit depth) a ZiX ({3, 7) were 
similar for all media tested. Furthermore, these stress 
profiles could be approximated (to the first order) by a 
simple scaling of generic stress profiles (supplementary 
text section 3, figs. S6 and S7, and tables S2 and S3). 

We next tested our hypothesis that forces on a com- 
plex intruder moving in granular media in the vertical 
plane could be approximated by the linear superposition 
of forces on all intruder elements. We measured the net 
lift F z and thrust F x on thin rigid model legs rotating 
about a fixed axle [simulating a tethered body (7, 11)] 
through granular media in the vertical plane at ~1 cm/s 
(Fig. 3 and movie S2) (28). We then compared them to 
predictions from the resistive force model by the integra- 
tion of stresses over the legs (movie S3) 
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where S is the leading surface of the leg; dA s1 \z\ s , (3 8 , 
and 7 S are the area, depth below the surface, angle of at- 
tack, and angle of intrusion of infinitesimal leg elements; 
and a ZjX ((3 s , j s ) are element stresses per unit depth (in- 
terpolated from data in Fig. 2, C and D). To test the 
robustness of our force model, we used three model legs 



of different geometries [with the same maximal leg length 
2R (28)]: a RHex robot's c-leg (21, 22), a flat leg, and 
a reversed c-leg (Fig. 3, A to C). In model calculations, 
each leg was divided into 30 elements. 

In all media tested (Fig. 3,D to F, and fig. S5), we 
observed that for all three legs, the measured net lift and 
thrust F ZyX as a function of leg angle (solid curves) were 
asymmetric to the vertical downward direction (0 = 0), 
and were larger during intrusion (0 < 0) than during ex- 
traction (0 > 0). Peak F ZjX were largest on the c-leg and 
smallest on the reversed c-leg. The reversed c-leg experi- 
enced significant negative lift (suction force, F z < 0) dur- 
ing extraction. For all media tested, our resistive force 
model predicted F ZjX versus for all three legs (dashed 
curves), capturing both the magnitudes and asymmetric 
profiles. The relative errors of peak forces between data 
and model predictions were within 10% for the c-leg in 
four media tested, and within 33% for all three legs in all 
media tested. The accuracy of our resistive force model 
was significantly better than that of previous force mod- 
els in which stresses depended only on depth (Fig. 3, D 
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FIG. 3. The resistive force model predicts forces on intrud- 
ers of complex morphology and kinematics moving in granular 
media (28). Three thin rigid model legs of different geome- 
tries, (A) a c-leg, (B) a flat leg, and (C) a reversed c-leg, 
were rotated about a fixed axle at a hip height h through 
granular media (gray area) in the vertical plane at an angular 
velocity cj, generating leg speeds of v ~ 1 cm/s, and net lift 
F z (blue) and thrust F x (green) were measured as a function 
of leg angle (movie S2). All three legs had identical maxi- 
mal length 2R from the axle, g is gravitational acceleration. 
(D to F) F z , x versus on the three legs measured in the ex- 
periment (solid curves) and predicted by the resistive force 
model (dashed curves) using Eq. 2 (movie S3). Insets: F z , x 
versus from experiment (solid curves) versus predicted (dot- 
ted curves) using Eq. 2 and previous force models in which 
stresses depended only on depth (supplementary text section 
1 and fig. S2). 
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to F, insets; supplementary text section 1 and fig. S2). 
Furthermore, our resistive force model revealed that the 
c-leg generated the largest forces, because its morphol- 
ogy allowed leg elements to not only reach deeper depths 
but also access larger stress regions in Fig. 2, C and D 
(particularly for elements at large depths). 

Our discovery of the insensitivity of the stress pro- 
files to particle properties (fig. S4) has practical ben- 
efits: For granular media of near-monodispersed, near- 
spherical, rounded particles, as an alternative to measur- 
ing a ZjX for all attack angles ft and intrusion angles 7 in 
the laboratory, one can simply perform a single measure- 
ment [of 0^(0,71-/2), using a horizontal plate penetrating 
vertically downward] to infer all a z ^ x (f3 : j) by a scaling 
routine (supplementary text section 3 and fig. S9) and 
predict forces (with a small loss in accuracy for the c- 
leg and the fiat leg, but a larger loss in accuracy for the 
reversed c-leg, fig. S8). 

We tested the ability of our resistive force model to 
predict legged locomotion. We chose to study the loco- 
motor performance (speed) of a small RHex-like robot 
(22) (Fig. 4A, top, and movie S4) moving on granular 
media (28). The robot's six legs rotated nearly entirely in 
the vertical plane during locomotion, and its small size 
ensured that leg intrusion speeds were low enough for 
particle inertia to be negligible. We chose poppy seeds as 
the test granular medium, because the grains were both 
small enough be prepared in our fluidized bed track (21) 
and large enough to not jam the robot's motor and gear 
trains. The robot's legs had a similar friction coefficient 
with poppy seeds to that of the model legs and were suffi- 
ciently rigid so that they experienced negligible bending 
during movement (28). 

Unlike the sand-swimming lizard, which moves within 
granular media quasistatically (thrust and drag are al- 
ways roughly balanced) (26), legged locomotion on the 
surface of granular media is dynamic (forces are not al- 
ways instantaneously balanced). As a result, the resistive 
force theory (which solves for speed by balancing forces) 
(6, 26) cannot be directly applied. Thus, to use our resis- 
tive force model to calculate robot speed, we developed a 
three-dimensional multibody dynamic simulation of the 
robot (Fig. 4A, bottom) (28). The simulated robot had 
the same body and leg morphology and used the same 
alternating tripod gait as the actual robot and had its 
motion constrained in the vertical plane. We divided 
each body plate and leg into 30 elements. The velocity 
v and angular velocity uj of the simulated robot's body 
were calculated by 



r v(t + dt)=v(t) + ^dt 
u(t + dt) uj(t) + fdt 



(3) 



medium, calculated from our resistive force model by the 
integration of stresses over each leg and the body using 
Eq. 2; m and / are the robot's mass and moment of 
inertia; and t and dt are time and time step. To test the 
robustness of our resistive force model and simulation, 
we used legs of seven geometries with different curvatures 
1/r (given maximal leg length 2R') (Fig. 4C, inset) and 
varied stride frequency / to up to 5 Hz (28). 

We observed similar robot kinematics (Fig. 4A) and 
forward speed v x versus time t (Fig. 4B) in both the 
experiment (movie S5) and simulation (movie S6). The 
robot moved faster and penetrated its legs less deeply 
during stance using c-legs (Fig. 4 A, left; Fig. 4B, red) 
than using reversed c-legs (Fig. 4A, right; Fig. 4B, blue). 
Average forward speed v x increased with stride frequency 
/ for legs of all curvatures 1/r and was lower at any / us- 
ing legs of negative curvatures than using legs of positive 
curvatures (Fig. 4, C and D). The agreement between 
experiment and simulation in v x (f, 1/r) was remarkable: 
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where F and N are the sum of net forces and torques 
on all the six legs and the body exerted by the granular 



FIG. 4. A multibody dynamic simulation using the resis- 
tive force model predicts legged locomotion on granular me- 
dia (28). (A) Side views of a small RHex-like robot (movie 
S4) at mid-stance during locomotion on granular media, us- 
ing c-legs (left) and reversed c-legs (right) in the experiment 
(top, movie S5) and simulation (bottom,movie S6). Arrows in 
the simulation indicate element forces on one tripod of legs. 
g is gravitational acceleration. (B) Forward speed v x versus 
time t from two representative runs using c-legs (red, stride 
frequency / = 2.0 Hz, curvature 1/r = 1/R') and reversed 
c-legs (blue, / = 2.2 Hz, 1/r = -l/R'). (C) Average forward 
speed v x versus / using legs of seven curvatures 1/r transi- 
tioning from reversed c-legs to c-legs (inset), where r is the 
radius of curvature, 2R' is the maximal length of the robot 
legs, and the minus sign denotes reversed legs. (D) v x versus 
1/r at / = 1, 2, 3, and 4 Hz. In (B) to (D), solid and dashed 
curves indicate the experiment and simulation, respectively. 
Error bars in (C) denote ±1 SD (<C 1 cm/s in the simulation). 
Experimental data in (D) are interpolated from those in (C) 
(hence no error bars). Red and blue arrows in (C) and (D) 
indicate averages from data shown in (B). 
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Errors were within 20% for 90% of the / and 1/r tested 
and within 35% for all the / and 1/r tested. Simulation 
using our scaling routine also achieved reasonable accu- 
racy (supplementary text section 5 and fig. SI 3). This 
was an improvement over simulation using previous force 
models in which stresses depended only on depth (fig. 
S13). Our resistive force model and simulation revealed 
that the robot moved faster using c-legs than using re- 
versed c-legs, because whereas the c-legs penetrated less 
deeply, their elements accessed larger stress regions in 
Fig. 2, C and D, resulting in larger leg lift (fig. S14) 
and smaller body drag. Our model and simulation also 
allowed the prediction of ground reaction forces on gran- 
ular media(Fig. 4A, red arrows, and fig. S13), which 
would be difficult to measure otherwise. Furthermore, 
our model and simulation predicted that using arc-like 
legs (given maximal leg length 2R r ) of an optimal curva- 
ture of 1/r = 0.86/i?', the robot would achieve maximal 
speed of v x = 72 cm/s 5 body length/s) at 5 Hz. 
Our approach affords significant reduction in the com- 
putational time needed to model movement on granular 
media. For example, relative to our multipart icle dis- 
crete element method (DEM) simulation of movement 
on granular media (23, 27), our simulation using the re- 
sistive force model can achieve a factor of 10 6 in speed-up 
(e.g., 10 s versus 30 days using DEM to simulate 1 s of 
locomotion on a granular bed of 5 x 10 6 poppy seeds). 

We close with a brief discussion of the limitations of 
our model. We tested the predictive power of our scaling 
routine (supplementary text section 3 and fig. S9) for 
two highly polydispersed, nonspherical, highly angular 
natural sands (supplementary text section 4, figs. S10 
to S12, and table S4). We found that the model accu- 
racy for natural sands was slightly worse than found in 
the glass spheres and poppy seeds (for example, 35% ver- 
sus 20% error in peak F z for a rotating c-leg). As was 
the case for the near-spherical granular media tested, the 
functional forms of forces on the c-leg and the flat leg 
were still well captured by our scaling routine. Further- 
more, the overestimation was not affected by reducing 
the polydispersity of the natural sand. This suggests that 
the nonsphericity and angularity of natural sand particles 
(30) may be the cause of this overestimation, which may 
require additional model fitting parameters and scaling 
factors. Our model is intended for dry sand [~0.1 to ~1 
mm in particle diameter (29)] and may not work for dry 
cohesive powder (<~0.01 mm) (31). We do not expect 
our model to work if the particle size approaches a char- 
acteristic length of the locomotor (for example, ~l-cm 
particles for our robot of ~l-cm foot size), so that the 
continuum assumption breaks down and particles effec- 
tively become boulders. We also do not expect our model 
to capture wet, cohesive flowable media such as soil and 
mud. 

We have developed a new approach to predicting 
legged locomotion on granular media. This terradynam- 



ics relies on new resistive force measurements and lin- 
ear superposition (6, 26, 27). The general profiles of 
these resistive force measurements are insensitive (other 
than magnitudes) to a variety of granular media com- 
posed of slightly polydispersed, approximately spherical, 
rounded particles. Our terradynamics may not be limited 
to legged locomotion, because the integration of stresses 
should in principle work for devices of other morphology 
and kinematics, such as wheels, tracks, and earth movers 
moving on granular media (2, 25). For the particle types 
tested here, an important addition to our model would be 
to capture three-dimensional effects (16) and spatial and 
temporal variation in compaction (21) and slopes (17), 
and test its validity in the high speed "inertial fluid" 
regime (when leg intrusion speed is > 1 m/s, at which 
particle inertia dominates forces) (23). Our resistive force 
model also provides opportunities to test and develop 
new physics theories of dense granular flow (32). Finally, 
we envision that, in concert with aero- and hydrodynam- 
ics (3-12), a general terradynamics of complex ground 
will not only advance understanding of how animals move 
(1) at present (5-10, 13-19, 26) and in the past (17, 33), 
but also facilitate the development of robots with loco- 
motor capabilities approaching those of organisms (11, 
12,20-23). 
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Materials and Methods 

Force measurements 

We used aluminum to construct the plate element (area 
A = 3.81 x 2.54 cm 2 , thickness = 0.64 cm) and model 
legs (maximal length 2R = 7.62 cm, width = 2.54 cm, 
thickness = 0.64 cm). We measured the friction coeffi- 
cient \i between aluminum and poppy seeds to be 0.40 
(later we constructed robot legs using plastic of a similar 
friction coefficient with poppy seeds, 0.36), by placing an 
aluminum plate on a wooden plate bonded with a single 
layer of poppy seeds, increasing the slope of the wooden 
plate from zero, and examining the angle £ at which the 
aluminum plate began to slide. Thus /i = tan£. The 
length and width of both the plate element and model 
legs were ~ 10 times the particle diameter, ensuring that 
the granular media could be approximated as a contin- 
uum. 

Before each force measurement, we used an air flu- 
idized bed (24 x 22 cm 2 surface area) (26) to prepare 
the granular media (15 cm deep) to a well-defined com- 
paction (see table [Si] for the volume fractions of closely 
and loosely packed states of the granular media tested). 
Air flow was turned off during force measurements. We 
used a 6 degree-of-freedom robotic arm (CRS Robotics) 
to move the plate element and rotate the model legs. 
We used a 6-axis force and torque transducer (ATI In- 
dustrial Automation) mounted between the intruder and 
the robotic arm to measure forces to a precision of 0.05 N 
at a sampling frequency of 100 Hz. We performed all the 
force measurements at low speeds (~ 1 cm/s) to ensure 
that particle inertia was negligible, and in a vertical plane 
at the middle of the air fluidized bed and far from the 
sidewalls (distance > 3 cm) to minimize boundary effects. 

In the plate element intrusion experiment, we attached 
the plate to the force and torque transducer via a sup- 
porting rod and an adjustable mount with which attack 
angle j3 could be varied. We varied intrusion angle 7 by 
adjusting the trajectory of the robotic arm. For each 
combination of f3 and 7, we separately measured the 
forces on the supporting rod and mount moving in the 
granular media without the plate, and subtracted them 
to obtain forces exerted by the granular media on the 
plate alone. 

During each test session, we first prepared the granular 
media while the plate was above the surface. We then 
moved the plate (oriented at attack angle /?) downward to 
the surface (depth \z\ = 0), paused it for 2 seconds, and 
then intruded it into the granular media along intrusion 
angle 7. After intrusion was finished, we prepared the 
granular media again, and extracted the plate along the 
same path. This gave us measurements of stresses a ZyX 
for both ±(/3, 7). For horizontal movements (7 = 0), 
&z,x were nearly constant when the plate was far from 
the container sidewalls, and we obtained a ZjX by fitting 



Eq. 1 to averages of a z ^ x in the steady state regions at 
three depths (\z\ = 2.54 cm, 5.08 cm, and 7.62 cm). We 
measured a ZyX for 7 within [— 7r/2, tt/2], and determined 
a ZjX for 7 within [— 7r, — tt/2] and [7r/2,7r] by symmetry: 

(a z (P,j)=a z (-P,-7r-j) if -7T < 7 < -tt/2 
\a z (P,'y)=a z (-P,7r-'y) if tt/2 < 7 < tt 

\a x (P,j) -OL x (-P, 7T 7) if -TT < 7 < -7T/2 

[a^(/3,7) = -a x (-P,TT- 7) if tt/ 2 < 7 < n 

(S2) 

In the model leg rotation experiment, we rotated the 
model legs at an angular velocity uj = 0.2 rad/s at a hip 
height h = 2 cm within leg angle — 37r/4 < < 37r/4 
(at the beginning and end of which all three legs tested 
were fully above the granular surface), where leg angle 
was defined as the angle sweeping from the vertical down- 
ward direction to the direction along which leg length was 
maximal. For each model leg, we separately measured 
the forces due to the weight of the model legs rotating 
in the air, and subtracted them to obtain the forces on 
the model legs exerted by the granular media during ro- 
tation. 

Due to the high repeatability of our fluidized bed 
and robotic arm, we found that for all media tested, 
run-to-run variation in a ZjX for fixed f3 and 7 was always 
within 0.005 N/cm 3 at any given depth; thus we only 
performed one trial for each combination of (3 and 7. 
All the stresses were calculated in the regions where the 
plate was far from the container boundaries (distance 
> 6 cm). We also confirmed that for low enough speeds, 
intrusion forces in granular media were insensitive to 
speed (for example, in loosely packed poppy seeds, at 
v = 1 m/s, force only increased by less than 20% from 
that at v = 1 cm/s); thus, particle inertia was negligible. 

Robot experiments 

We built our robot (body length = 13 cm, body mass 
= 150 g) by modifying a small commercially available 
robot (RoboXplorer, Smart Lab). The robot had simi- 
lar morphology and kinematics as a RHex robot, with a 
rigid body and six legs performing 1 degree-of-freedom 
rotation in an alternating tripod gait. We substituted 
the stock motor with a stronger one (RadioShack Super 
Speed 9-18 VDC Hobby Motor, Model # 273-256) and 
modified the gear trains (gear ratio: 47 revolutions in 
the motor transmits into 1 rotation of the legs). These 
changes increased maximal stride frequency / to 5 Hz. 
We removed the external body shell to reduce weight 
and the belly area (to 13 x 2 cm 2 ). This reduced drag on 
the belly during locomotion on granular media. 

We used 3-D printing to make custom robot legs. All 
the legs had the same maximal length 2R! = 4.1 cm, 
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width = 1.0 cm, and thickness = 0.3 cm, but different 
curvatures 1/r = [-1, -0.92, -0.60, 0, 0.60, 0.92, l]/i?'. 
The ABS plastic used to fabricate the legs had a sim- 
ilar friction coefficient with poppy seeds (0.36) to that of 
aluminum with poppy seeds (0.40). The leg width and 
length were ~ 10 times the particle diameter, allowing 
the granular media to be approximated as a continuum. 
We ensured that the legs had large enough stiffness and 
moved like rigid bodies (< 5% deformation) during loco- 
motion. 

We tuned the center of mass of the robot to overlap 
with the geometric center of the body by adding mass to 
the lighter end of the robot. We measured the masses, 
dimensions, and relative positions of all robot body and 
leg parts, and calculated the moment of inertia of the 
robot to be / = 2.08 x 10 3 g cm 2 about the pitch axis 
through the center of mass. We powered the robot by 
an external power supply (Power Ten Inc.) to ensure 
constant voltage during trials, and adjusted voltage to 
vary stride frequency / between trials. 

Before each trial, we used an air fluidized bed track 
(200 x 50 cm 2 surface area) to prepare the granular 
media (12 cm deep) to a well-defined compaction, using 
methods similar to those in (21). We used two synchro- 
nized high-speed cameras (X-PRI, AOS Technologies) 
to capture top and side views of the robot's locomotion 
at 500 frame/s. We measured / from the side view, 
and measured forward speed v x from the top view by 
digitizing a high contrast marker placed near the center 
of mass. For each / and 1/r, we performed three trials 
and reported mean ± s.d. for average forward speed v x 
in experiment. 



stride frequency / between < / < 5 Hz in increments 
of 0.2 Hz. For each / and 1/r, we performed three trials 
using different initial conditions with a phase difference 
of 27r/3 for each tripod. We found that this resulted in 
variation in average forward speed v x of <C 1 cm/s; thus, 
we reported only the means of v x in simulation. We con- 
firmed that at the maximal stride frequency tested (5 Hz) 
the leg speeds were < 1 m/s averaged over a stance, al- 
lowing particle inertia to be negligible. 



Robot simulation 

We used a multibody dynamic simulator, MBDyn (34), 
to create a simulation of the robot locomotion on granu- 
lar media. MBDyn features a full three-dimensional sim- 
ulation with 6 degrees of freedom (3 translations and 3 
rotations). We constructed the robot body using 3 rigid 
plates (the front, rear, and belly surfaces) and 6 rigid 
legs. We constrained the robot body movement within 
the vertical plane, with fore-aft and dorso- ventral transla- 
tions and pitch, and allowed the legs to only rotate about 
their axles perpendicular to the vertical plane. For each 
time step, we calculated the depth \z\ s , attack angle /3 S , 
and intrusion angle 7 S for each element to determine el- 
ement stresses a Z}X (using Eq. 1 and a Z}X interpolated 
from the data in Fig. 2, C and D). We then summed 
forces on all elements to obtain net forces F ZjX using our 
resistive force model by Eq. 2, and calculated the body 
dynamics by Eq. 3. We found that dividing each body 
plate and leg into fewer elements could further increase 
simulation speed at the cost of model accuracy. 

In simulation tests, we varied curvature 1/r between 
-l/R' < 1/r < l/R' in increments of 0.04/i?', and varied 
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Supplementary Text 

This supplementary text contains 5 sections: 

In section 1, we review the depth-only vertical pene- 
tration and horizontal drag force models. Based on the 
flat-plate approximation, we then calculate stresses per 
unit depth and net forces on rotating legs using these 
depth-only force models, and compare with our resistive 
force measurements and model predictions. 

In section 2, we present our resistive force measure- 
ments and model calculations for a variety of granular 
media tested, which have different particle size, density, 
friction, and compaction; from these results we discover 
that the stress profiles are generic. 

In section 3, we develop a scaling routine to capture 
the similar stress profiles observed for the variety of gran- 
ular media, by scaling generic stress profiles determined 
from averages of fits to the measured stress profiles. This 
provides a practical means to easily apply our resistive 
force model. 

In section 4, we test the ability of our scaling routine 
to predict forces on intruders moving in natural sands, 
and discuss limitations of our model for natural sands. 

In section 5, we compare the predictive accuracy for 
locomotor performance (speed) of the robot using our 
resistive force model to that using our scaling routine 
and that using the depth-only force models. 

1. Depth-only force models and their limitations 

Previous studies of intrusion forces in granular me- 
dia focused on simple intruders (e.g. plates, rods, and 
spheres) moving with simple kinematics (e.g. vertical 
penetration and horizontal drag). For example, the verti- 
cal force f z on a horizontal plate element (f3 = 0) moving 
vertically (7 = ±7r/2) in granular media was observed to 
be proportional to the plate's depth \z\ and area A (21, 
35): 



f z = a z (Q,8ga(z)ir/2)\z\A 



(S3) 



Similarly, the horizontal force f x on a vertical plate ele- 
ment (/? = ±7r/2) moving horizontally (7 = 0) in granu- 
lar media was proportional to plate depth \z\ and plate 
area A (36-38): 



f x =a x (ir/2,0)\z\A 



(84) 



where z is the velocity of the plate in the vertical di- 
rection, and a z (0, sgn(i)7r/2) and a^(7r/2,0) are verti- 
cal and horizontal stresses per unit depth for (/?, 7) = 
(0, ±7r/2) and (7r/2,0) (determined from measurements 
in Fig. 2, C and D). 

Both the vertical penetration and horizontal drag force 
models only account for the dependence of stresses on 
the intruder's depth (|z|), but not the dependence on 
its orientation (attack angle /?) or movement direction 




afl (N/cm 3 ) a x 



jo, ^yu^gi 




FIG. SI. Effective vertical (a[ z ^(/3,j), left) and horizon- 
tal (a^i(/?,7), right) stresses per unit depth as a function 
of attack angle f3 and intrusion angle 7 for loosely packed 
poppy seeds, calculated from eq. |S7| using a z (0, sgn(z)n /2) 
and 0^(71-/2,0) measured in experiment (Fig. 2, C and D). 
See Fig. 2, A and B for schematic of the experiment and def- 
inition of variables. 



C-leg 



Flat leg 



Reversed c leg 



Poppy seeds 
loosely packed 



Poppy seeds § 
closely packed y£ 



0.3 mm glass spheres 
loosely packed LL 



0.3 mm glass spheres 
closely packed 



3 mm glass spheres * 
closely packed ^ 




M 



N 











/ 
/ 

J< 


\ 
\ 

\ 







*J2 -n/2 



n/2 -n/2 



FIG. S2. Net lift F z (blue) and thrust F x (green) versus leg 
angle on the three model legs for all media tested. Solid 
curves: experimental measurements. Dotted curves: predic- 
tions from the depth-only force models. See Fig. 3 A to C for 
schematic of the experiment and definition of variables. 



(intrusion angle 7). Hereafter we refer to these force 
models as the "depth-only force models" . 

In previous studies of legged locomotion on granular 
media (16, 21, 37), due to the lack of the resistive force 
model, the forces on a complex intruder were estimated 
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from the depth-only force models, using the flat-plate 
approximation: The vertical force on a leg element of area 
dA s and arbitrary f3 and 7 was approximated by that on 
a horizontal leg element whose area was the element area 
projected into the horizontal plane (16, 21): 



model legs rotated through all media tested, by the inte- 
gration of stresses over the legs (Eq. 3) using ai*i(/?,7). 
We found that they not only significantly overpredicted 
peak F z (by up to 1400%) but also erroneously predicted 
similar peak F zx for the c-leg and reversed c-leg (fig. S2). 



dA zs = | cos/3 |cL4 s 



(S5) 



Similarly, the horizontal force on the leg element was ap- 
proximated by that on a vertical leg element whose area 
was the element area projected into the vertical plane 
(37): 



2. Generality of resistive force model for a variety of 
granular media 

To test the generality of our resistive force model and 
provide a database for future studies (we provide the 
force measurements in this section in Additional Data 



dA zs = |sin/3|cL4 s 



(S6) 



The effective stresses per unit depth calculated from 
these depth-only force models using flat plate approxi- 
mation were then: 



Mi 



az\Pn) = a*(0,sgn(i)7r/2)|cos/3| 
,7) = ^(7r/2,0)|sin/3| 

M 



Poppy seeds 
loosely packed 



(S7) 



Comparing the ai,]c(/3,7) calculated using eq. 



S7 



(fig. SI) to our resistive force measurements a ZyX (/3^) 



Poppy seeds 
closely packed 



(Fig. 2, C and D) for loosely packed poppy seeds, we 
found that these depth-only force models did not cap- 
ture most of the measured stress profiles. In particular, 
the depth-only vertical penetration force model overpre- 
dicted a z (/3, 7) for all /3 and 7 except when the plate was 
horizontal (/3 = 0) and moving vertically (7 = ±7r/2). 

\z\ 

Contrary to experimental measurements, the a Z)X cal- 
culated from the depth-only models were symmetric to 

M 



P = (horizontal orientation), and a l z l (or a l x l ) was al- °"^etpacktT P a ^ 



ways opposing the plate's vertical (or horizontal) velocity 
(the shaded regions shrank to a line). 

We used these depth-only force models to calculate the 
net lift F z and thrust F x versus leg angle 6 for the three 
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FIG. S3. The three granular media tested in our study. Top: 
regular images. Bottom: microscope images. The length of 
each scale bar is 1 mm. Photo credit of (A) and (C): Sarah 
Sharpe. 



FIG. S4. Vertical (a[ zl (/3, 7), left) and horizontal (al*i(/3, 7), 
right) stresses per unit depth as a function of attack angle f3 
and intrusion angle 7 for all media tested. (A) and (B) are 
reproduced from Fig. 2, C and D. See Fig. 2, A and B for 
schematic of the experiment and definition of variables. 
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FIG. S5. Net lift F z (blue) and thrust F x (green) versus leg 
angle on the three model legs for all media tested. Solid 
curves: experimental measurements. Dashed curves: resis- 
tive force model predictions. (A) to (C) are reproduced from 
Fig. 3, D to F. See Fig. 3 A to C for schematic of the experi- 
ment and definition of variables. 



Table S5), we performed resistive force measurements 
for three granular media — poppy seeds, 0.3 mm glass 
spheres, and 3 m glass spheres (fig. S3), which have 



different particle size, shape, density, and friction (mea- 
sured by angle of repose) (table [Si]). We prepared these 
granular media into well-defined compactions which af- 
fected stresses (21, 38). We prepared poppy seeds and 
0.3 mm glass spheres into both a loosely packed (LP) 
and a closely packed (CP) states, and prepared 3 m glass 
spheres into a closely packed (CP) state (table [Si]). 

Despite differences in magnitudes and fine features, we 
observed similar profiles of the vertical and horizontal 
stresses per unit depth a ZyX {P^) for all media tested 
(fig. S4). We provide the a z ^ x (f3,j) data for all media 



in Additional Data Table S5 (separate file in Microsoft 
Excel format). 

For all media tested, our resistive force model pre- 
dicted the net lift F z and thrust F x versus leg angle 6 
on the three model legs rotated through granular media 
(fig. |S5[). Compared with predictions from the depth- 



only force models (fig. SI), our resistive force model had 
a significant improvement in accuracy. 

3. Scaling routine for easy use of the resistive force 



model 

To provide a means for practical use of our resistive 
force model and comparison with new theories of dense 
granular flow (32), we performed a fitting approximation 
to the stress per unit depth data a ZyX for all media tested, 
and developed a scaling routine based on the data fits. 

We first performed a discrete Fourier transform of the 
<*z,x(P,l) data (fig. [Si]) over -tt/2 < /3 < ir/2 and -tt < 



7 < 7r to obtain a fitting function. We examined the 
Fourier coefficients and found that the a z ^ x (f3,j) data 
of all media tested could be well approximated by the 
following fits: 



a 



fit 



.fit 



08,7)= Yl ^2[^m jn cos27r( 

m= — l n=0 

+B rn , n sin27r( 
l l 

(&7)= E Y,[Cm,nCos27r( 

=0 

-\-D rn ^sin2ir( 



-1 n=0 



m/3 

TT 

mf3 

TT 

mf3 

TT 
TT 



(S8) 



717 
~ 2tt 

2ir J 



using nine zeroth- and first-order terms (whose magni- 



tudes are larger than 0.05A ,o) (table S2): 



M = (^o,o ^i,o #1,1 #o,i #-i,i Ci,i Co,i C-i,i #i,o) 

(S9) 

By symmetry, a z (/3 < 0,7r/2) = a z (f3 > 0,7r/2), and 
ot x (fi < 0,7r/2) = — a x (P > 0,7r/2). However, the data 
slightly deviated from this equality because the initial 
positions of the plate were close to one of the boundaries 
of the container. Therefore, before the Fourier transform, 
we averaged the raw data for 7 = ±tt/2 to restore sym- 
metry by a z (fi < 0,7r/2) = a z (/3 > 0,tt/2) = \[a z (l3 < 
0, 7r/2) + a z (f3 > 0,7r/2)]. The raw data were also not 
uniformly sampled in the 7 direction; we found that this 
only resulted in small errors in data fitting. 
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FIG. S6. Scaled Fourier coefficients Mi/Q of all media tested 
(thin colored curves) can be approximated by a generic coef- 
ficient curve Mo (thick black curve). 
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FIG. S7. Generic stress (per unit depth) profiles 

af?r ric (A7) for all media tested. See Fig. 2, A and B for 
schematic of the experiment and definition of variables. 



We found that for all media tested (denoted by i = 
poppyLP, poppyCP, 0.3mmLP, 0.3mmCP, and 3mmCP), 
the Fourier coefficients Mi could be collapsed onto a 
generic coefficient curve, Mo, by dividing Mi by a scaling 



factor d (table [S2J fig.|S10j): 

Mi/Ci 



M 



(S10) 



This enabled us to scale stresses per unit depth (a ZjX ) 
and thus forces (f ZyX and F ZyX ) for all media tested. By 
eq.|S8| using the generic coefficient curve Mq (£ = 1) from 



table |S2[ we calculated generic stress (per unit depth) 
profiles af£ eric (/3, 7 ) (fig.|S7|). We found that, by multi- 



plication by the scaling factor £, these generic stress pro- 
files well approximated the measured a ZjX (l3,j) (fig. S4) 
for all media tested: 



c<r ric (/^,7) 



(Sll) 



Using Eq. 2 and a|^ eric (/3, 7), we calculated the 
generic force profiles Ff^ neric (#) on the three model legs 
rotated through granular media (fig. |S8[). We found that 
in all media tested and for both the c-leg and the fiat 
leg, these generic force profiles captured the measured 
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F z ^ x scaled by the scaling factors £: 

« F z , x (0)/C 



(S12) 



However, the agreement was worse for the reversed c-leg. 
Our resistive force model revealed that this was because 
by using Ff^ neric (#), stresses were significantly overesti- 
mated for the reversed c-leg's elements that reached large 
depths. 

We further observed that the ratio of the maximal ver- 
tical stress (which occurred on a horizontal plate mov- 
ing downward) between the measurements (a z (0,7r/2)) 
and fits (0^(0, tt/2)) was similar for all media tested (ta- 
ble [Si: 



X = a z (0, 7r/2)/af (0, tt/2) = 1.26 ± 0.14 (S13) 

where x is i n mean ± s.d. 

Therefore, we propose that for a sufficiently level 
and uniform dry granular medium composed of near- 
monodispersed, near-spherical, rounded particles of ~ 
0.1 to ~ 1 mm in diameter, one can simply measure its 
maximal vertical stress a 2 (0,7r/2) by pushing a horizon- 
tal plate downward to infer the maximal value of the fit 
vertical stress a^ t (0, 7r/2): 

af (0, tt/2) = a z (0, tt/2)/x « 0.8a*(0, tt/2) (S14) 

The value of 0^(0, tt/2) in N/m 3 then gives the scaling 
factor for this granular medium, because we choose the 
magnitudes of the generic curve Mo so that the values 
of ( and 0^(0,71-/2) are nearly the same for all media 
tested: 



a 



fit 



08. 7) 



(S15) 



Note that this equation only equates the numeric values 
on both sides, because £ is dimensionless. 

Then, from eq . |S11[ by scaling the generic stress pro- 
files «| e " eric (fig. S7) by the determined scaling factor 



Push a horizontal plate downward to measure tx z (0, tt/2) 



0^(0, n/2) * 0.8 a z (0, tt/2) 



Scaling factor £ « a^ x {0, n/2) 



a ZA (ft Y) - a» w (ft y) = ^ af e x neric (p, y) 



FIG. S8. The measured net forces scaled by the mea- 
sured scaling factors F Z , X (Q)/^ for all media tested (thin 
curves), in comparison with the generic profiles of net forces 
Fj^ eric (0) (thick curves) calculated from the generic stress 
profiles af^ eric (%-[S7l)- See Fig. 3 A to C for schematic of 
the experiment and definition of variables. 



Use a z c x aled (p, y) to calculate forces and dynamics by Eqs. 2 and 3 



FIG. S9. A scaling routine to easily apply our resistive force 
model. 
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we can obtain an approximation of stress profiles for this 
granular medium: 



29 Palms sand 



29 Palms sand 
(reduced polydispersity) 



a. 



,x(/3,7) 



a' 



scaled 



(/? )7 )=C<r ric (/?,7) (S16) 



This scaling routine provides an alternative to measur- 
ing a ZjX for all /3 and 7 (at the cost of model accuracy). 
As demonstrated by the model leg rotation experiments 



(fig. S8 ) and robot locomotion experiments (see fig. |S13| 
in the next section), our scaling routine only suffers a 
small loss in accuracy for much of the leg morphology 
and stride frequencies tested. This technique can be par- 
ticulary useful in a field setting, because only a single 
force measurement is needed. 

We summarize these practical steps to use our resistive 
force model in fig. [S9 



4. Applicability of resistive force model to natural sands 

To test the predictive power of our resistive force model 
for natural sands, we chose two natural sands of higher 
polidispersity and angularity than those of the granular 
media tested, and examined whether the scaling routine 
(fig. |S9|) could predict the net forces F Z:X (6) on the three 
model legs during rotation through the natural sands. 

The two natural sands, referred to as "Yuma sand" and 



"29 Palms sand" (table |S4} fig. S10), were collected from 
the Mojave Desert in the southwest United States, one 
from Yuma, Arizona and the other from 29 Palms, Cali- 
fornia. The Yuma sand had most particles (70% by mass) 
in the 0.06-0.3 mm particle size range, and the 29 Palms 
sand had most particles (91% by mass) in the 0.3-3 mm 
particle size range. Both naturals sands were fluidized by 
the fluidized bed before each force measurement was per- 
formed. We further tested modified 29 Palms sand with 
reduced polydispersity to examine the effect of polydis- 
persity on model accuracy. 



29 Palms sand 



29 Palms sand 
(reduced polydispersity) 

cT 




FIG. S10. The natural sands used to test the predictive 
power of our scaling routine. Top: regular images. Bottom: 
microscope images. The length of each scale bar is 1 mm. 



A / 




B / 




c 






5 




5 





















|z| (cm) 



FIG. Sll. Maximal vertical stress <j z (0,tv/2) versus depth 
\z\ measured in the natural sands using a plate element 
horizontally oriented and penetrating vertically downward 
(/3 = 0, 7 = 7r/2). Dashed lines are linear fits with zero inter- 
cept to the data in the linear regime at shallow depths. See 
Fig. 2, A and B for schematic of the experiment and definition 
of variables. 



During downward penetration (7 = tt/2) into both the 
Yuma and 29 Palms sands, the vertical stress <j z on a 
horizontal plate (/? = 0) increased nearly proportion- 
ally to depth \z\ at low enough depths (fig. Sll, solid 
curves), similar to observations in the granular media 
tested. However, for both the natural sands, <j z vesus \z\ 
displayed nonlinearity as depth further increased sooner 
than observed for the granular media tested. This was 
because both natural sands had larger internal stresses 
(larger a ZyX magnitudes) and likely suffered boundary ef- 
fects at much shallower depths as compared to the gran- 
ular media tested. In addition, for the Yuma sand which 
had more particles in the 0.06-0.3 mm particle size range, 
this nonlinearity due to boundary effects was more pro- 
nounced. Therefore, to minimize possible errors from 
boundary effects, we perform linear fits with zero inter- 



C-leg 



Flat leg 



Reversed c-leg 



c 

-Exp 

-- Scaling (u 


sing Cag neric ) 







29 Palms sand 
(reduced polydispersity) n 10 




-ti/2 ti/2 -Tt/2 Tt/2 -jt/2 Tt/2 



FIG. S12. Net lift F z (blue) and thrust F x (green) ver- 
sus leg angle on the three model legs for the natural sands 
tested. Solid curves: experimental measurements. Dashed 
curves: resistive force model predictions using the scaling rou- 
tine (fig. |S9| ). See Fig. 3 A to C for schematic of the experi- 
ment and definition of variables. 
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cept to the a z versus \z\ data (fig. |S11| dashed curves) 
over the linear regime (0-2 cm for the Yuma sand; 0- 
3 cm for the 29 Palms sand) to obtain a z (0, tt/2) (slopes 
of dashed lines) for both the natural sands. This gave 
us their scaling factors £ and estimates of th eir s tress 
profiles from a^ led (/3, 7) = (af£ eric (f3,>y) (eq.[S16|. 



For the three model legs rotated through both the 
Yuma and 29 Palms sands, both the net lift F z and 
thrust F x versus leg angle displayed asymmetric profiles 
similar to those observed for the granular media tested 
(fig. |S12 , solid curves). Using the a 2 (0,7r/2) obtained 
from the vertical penetration experiment (fig. |S11| ), our 
scaling routine significantly overpredicted F ZyX (6) on the 
c-leg for both the Yuma and 29 Palms sands (both by 
35%). Nevertheless, the shape of model predictions were 
similar to experimental observations (fig. |S12| dashed 



curves). The model accuracy was best for the fiat leg, 
and worst for the reversed c-leg. 

There are two major differences between the natural 
sands and the granular media tested: The natural sands 
have higher polydispersity, and are also less spherical and 
more angular. To probe where the overprediction (of 
forces on the c-leg) stemmed from, we further tested the 
29 Palms sand with reduced polydispersity (obtained by 



Resistive force model 



Scaling routine 



Depth-only models ^ 




sieving the sand to obtain only the 0.6-0.7 mm particles). 
We found that while both a z (0, tt/2) and F z ^ x (6) dropped 
by 14% (figs. Sll and S12), the model overprediction of 
forces on the c-leg remained the same (35%). This sug- 
gested that the non-spherical and angular shape of the 
natural sand particles, rather than their higher polydis- 
persity, was likely the cause of the model overprediction. 

5. Comparison of model accuracy in predicting legged 
locomotion 

To evaluate the accuracy of our scaling routine in pre- 
dicting legged locomotion, we used our multibody dy- 
namic simulation to simulate the robot's movement on 
loosely packed poppy seeds using the scaled generic stress 
profiles (af* led = (af™ eric ). We calculated the robot's 
average forward speed v x versus stride frequency / and 



leg curvature 1/r (fig. S13), and compare it with exper- 
imental measurements and resistive force model predic- 
tions (without using the scaling routine). As a compar- 
ison, we also calculated v x (f,l/r) using a^ z } x obtained 
from the depth-only force models (fig. SI). 

We found that the simulation using scaled generic 
stress profiles (C a f e £ enc ) on ly suffered a small loss in 
accuracy as compared to that using the measured a ZjX . 

\z\ 

By contrast, the simulation using stresses or z , x from the 
depth-only force models suffered a larger loss in accuracy 
(by up to 150% error), and in particular, erroneously pre- 
dicted that the robot moved at similar speeds using both 
c-legs and reversed c-legs. 

In addition to the predictions locomotor kinematics, 
our resistive force model and multibody dynamic simu- 



A 4 



C-leg 

Reversed / 
c-leg / 


* Tripod 1 n 

\ Tripod 2 » » 

I 1 1 
\ 1 1 
\ 11 





f(Hz) 




FIG. S13. Average forward speed v x of the robot as a func- 
tion of stride frequency / (left) and leg curvature 1/r (right) 
on loosely packed poppy seeds. Solid curves: experimental 
measurements. Dashed curves: simulation predictions. Top: 
using the measured stresses (a z , x )- Middle: using the scaled 
generic stress profiles (a s z c ^ led = CaT£ eric ) . Bottom: using 
stresses from the depth-only force models (ai*i). The drop 
of v x vs. / at near R' /r — 1 using ai*i was due to pitch in- 
stability of the robot in simulation (the robot became upside 
down). (A) and (B) are reproduced from Fig. 4, C and D. 
See Fig. 4 A and B for representative runs of the experiment 
and definition of variables. 



1 1 ■ 1 

0.75 1 1.25 

f(s) 



FIG. S14. Vertical (F z , top) and horizontal (F x , bottom) 
ground reaction forces versus time t on both tripods of the 
robot during locomotion on loosely packed poppy seeds, cal- 
culated from the simulation using the resistive force model. 
Solid and dashed curves correspond to forces on the two 
tripods of legs, respectively. Results are for the two repre- 
sentative runs shown in Fig. 4B. Note that these ground re- 
action forces do not include those on the robot body (which 
were also calculated but not plotted here). 
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lation also allowed us to calculate ground reaction forces 
(fig. |S14[ ) during locomotion on granular media (which 
would be otherwise difficult to measure), in a similar 
fashion as calculations of lift, drag, thrust, traction, and 
resistance in studies of vehicle mobility in fluids and ter- 
ramechanics studies. Notice the larger vertical ground 
reaction forces produced by the C-legs relative to those 
by the reversed c-legs. This allowed the robot to main- 
tain higher lift above the granular medium, reducing drag 
on the belly. 
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TABLE SI. Physical properties of the three granular media (in different compactions) for which resistive forces were measured. 
*Angle of repose for 3 mm glass spheres courtesy of Sarah Sharpe. 



granular medium 


particle diameter (mm) 


particle material density (g/cm 3 ) 


compaction 


volume fraction 


angle of repose (°) 


poppy seeds 


0.7 ±0.2 


1.1 


loosely packed (LP) 
closely packed (CP) 


0.58 
0.62 


36 
47 


0.3 mm glass spheres 


0.27 ±0.04 


2.5 


loosely packed (LP) 
closely packed (CP) 


0.58 
0.62 


25 
35 


3 mm glass spheres 


3.2 ±0.2 


2.6 


closely packed (CP) 


0.63 


21* 



TABLE S2. Zeroth- and first-order Fourier coefficients M (in N/cm 3 ) for all media tested. These coefficients are > 0.05Ao,o- 
Mo is a generic coefficient curve to which the M for each granular medium can be collapsed onto by division of a scaling factor 
C (C — 1 f° r Mo). We choose the magnitudes of M such that for all media tested the values of ( are nearly the same as the 
values of af t (0,7r/2) (this becomes useful later in the scaling routine). CP and LP indicates closely packed and loosely packed 
states. 



granular medium 


poppy seeds 


0.3 mm glass spheres 


3 mm glass spheres 


generic coefficients 


compaction 


LP 


CP 


LP 


CP 


CP 


n/a 


matrix notation 


MpoppyLP 


M p0 ppyCP 


Mo.3mmLP 


Mo.3mmCP 


M.3 mm CP 


Mo 


A ,o 


0.051 


0.094 


0.040 


0.081 


0.045 


0.206 


v4i, 


0.047 


0.092 


0.030 


0.063 


0.031 


0.169 




0.053 


0.092 


0.045 


0.078 


0.046 


0.212 




0.083 


0.151 


0.077 


0.133 


0.084 


0.358 




0.020 


0.035 





0.024 


0.012 


0.055 


Ci,i 


-0.026 


-0.039 


-0.031 


-0.050 


-0.031 


-0.124 


Co,i 


0.057 


0.086 


0.066 


0.087 


0.060 


0.253 


C-i,i 





0.018 











0.007 


#1,0 


0.025 


0.046 


0.017 


0.033 


0.015 


0.088 


scaling factor £ 


0.248 


0.488 


0.194 


0.371 


0.214 


1 



TABLE S3. Comparison of the measurements and fits of maximal vertical stress per unit depth (in N/cm 3 ) for all media 
tested. CP and LP indicates closely packed and loosely packed states. 



granular medium 


poppy seeds 


0.3 mm glass spheres 


3 mm glass spheres 


compaction 


LP 


CP 


LP 


CP 


CP 


a*(0,7r/2) 


0.35 


0.56 


0.24 


0.40 


0.29 


af (0,tt/2) 


0.26 


0.47 


0.19 


0.38 


0.22 


X = a,(0,7r/2)/a^(0,7r/2) 


1.37 


1.19 


1.27 


1.05 


1.33 



TABLE S4. Particle size distribution of the natural sands tested. 



natural sand 


particle diameter (mm) 


mass percentage (%) 




< 0.06 


2 


Yuma sand 


0.06-0.3 


68 




0.3-3.0 


17 




>3.0 


13 




< 0.3 


5 




0.3-0.6 


36 


29 Palms sand 


0.6-0.7 


55 




0.7-3.0 


3 




>3.0 


1 


29 Palms sand (reduced polydispercity) 


0.6-0.7 


100 
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Movie SI 

http : //crablab .gatech. edu/pages/movies/f iles/ 
Terrady n_sl .mov| 

Plate element intrusion experiments. Left: Intru- 
sion for representative attack angle and intrusion angle 
(/?, 7) = (7r/6,7r/4). Right: Extraction for representative 
(/?, 7) = — (tt/6, 7r/4). Top: Videos of the plate mov- 
ing in a granular medium (loosely packed poppy seeds). 
Bottom: Vertical (a z , blue curve) and horizontal (<r x , 
green curve) stresses versus depth \z\. Videos were taken 
with the plate at the container boundary for illustration 
purpose only. Forces were measured with the plate in 
the middle of the container, far away from the boundary. 
The red lines originating from the center of the top panels 
indicate granular forces on the plate element. The plate 
was moved at 1 cm/s. See Fig. 2, A and B for schematic 
of the experiment and definition of variables. 

Movie S2 



http : / / crablab . gatech . edu/ pages/movies/ f iles/ 
Terradyn _s2 . mov| 

Leg rotation experiments. Left: C-leg. Right: Re- 
versed c-leg. Top: Videos of the leg rotating through a 
granular medium (loosely packed poppy seeds). Bottom: 
Net lift (a z , blue curve) and thrust (cr Xl green curve) ver- 
sus leg angle 0. Videos were taken with the plate at the 
container boundary for illustration purpose only. Forces 
were measured with the plate in the middle of the con- 
tainer, far away from the boundary. The red lines origi- 
nating from the center of the top panels indicate granular 
forces on the model legs. The model legs were rotated 
at 0.2 rad/s (leg tip speed ~ 1 cm/s). Note that in the 
c-leg video the axle deflected slightly during rotation — 
this was an artifact due to the extra-long supporting rod 
in order for the plate to reach the boundary for illus- 
tration purpose only. For force measurements performed 
in the middle of the container, the supporting rod was 
shorter and the axle did not deflect. See Fig. 3 A to C 
for schematic of the experiment and definition of vari- 
ables. 

Movie S3 



http : / / crablab . gatech . edu/ pages/movies/ f iles/ 
Terradyn_s3 . mov 



Leg rotation model calculations. Left: C-leg. Right: 
Reversed c-leg. Top: Videos of the leg rotating through a 
granular medium (loosely packed poppy seeds). Bottom: 



Net lift (cr z , blue curve) and thrust (a x , green curve) 
versus leg angle 6. The thick red lines originating from 
the center of the top panels indicate net forces on the 
model legs. The thin red lines on the model legs indicate 
element forces (on a larger scale). See Fig. 3 A to C for 
schematic of the experiment and definition of variables. 
Movie S4 



http : //crablab . gatech . edu/pages/movies/f iles/ 
Terradyn _s4 . mov| 

Robot (body length = 13 cm, body mass m = 150 
g) running on a granular medium (loosely packed poppy 
seeds) using c-legs and stride frequency / = 5 Hz. Video 
is played in real time. 

Movie S5 



http : / / crablab . gatech . edu/ pages/movies/ f iles/ 
Terradyn_s5 .mov 

Robot experiments. Top: Using c-legs (/ = 2.0 Hz). 
Middle: Using reversed c-legs (/ = 2.2 Hz). Bottom: 
Forward speed v x versus time t. Videos are played 5 
times slower than real time. Note that in the trial using 
reversed c-legs, the robot body orientation was reversed; 
we verified that this did not affect robot kinematics and 
speed (because the center of mass of the robot was tuned 
to overlap with the geometric center of the body). See 
Fig. 4 A for definition of variables. 

Movie S6 



http : / / crablab . gatech . edu/ pages/movies/ f iles/ 
Terradyn _s6 .mov| 

Robot simulation. Top: Using c-legs (/ = 2.0 Hz). 
Middle: Using reversed c-legs (/ = 2.2 Hz). Bottom: 
Forward speed v x versus time t. Videos are played 5 times 
slower than real time. The red and green arrows on the 
two tripod of legs indicate the ground reaction forces on 
the leg elements calculated from simulation. Note that 
the forces on the body are not shown. See Fig. 4 A for 
definition of variables. 

Additional Data Table S5 (separate file in 
Microsoft Excel format) 



http : //crablab . gatech . edu/pages/movies/f iles/ 
Terrady n_table_s5 . xls| 



Original data of vertical (a z ) and horizontal (a x ) stresses 
per unit depth versus attack angle j3 and intrusion angle 
7, for all media tested. These data were used to produce 



Fig. 2, C and D and fig. S4 



